options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]' 
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -412584 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       37       -38.594    0.00664654    0.00175984           1           1       93    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -38.59
##    b0         5.57
##    b1         2.05
##    s          4.18
##    m0[1]      7.76
##    m0[2]     28.44
##    m0[3]     29.48
##    m0[4]     11.26
##    m0[5]     33.90
##    m0[6]     21.13
##    m0[7]     20.55
##    m0[8]     42.69
##    m0[9]     17.27
##    m0[10]    30.74
##    m0[11]    17.52
##    m0[12]    18.34
##    m0[13]    31.84
##    m0[14]    17.17
##    m0[15]    30.77
##    m0[16]    19.29
##    m0[17]    30.10
##    m0[18]    17.60
##    m0[19]    12.71
##    m0[20]    32.29
##    y0[1]      2.17
##    y0[2]     27.18
##    y0[3]     29.98
##    y0[4]     10.35
##    y0[5]     31.45
##    y0[6]     19.69
##    y0[7]     21.51
##    y0[8]     41.30
##    y0[9]     21.35
##    y0[10]    29.66
##    y0[11]    12.77
##    y0[12]    15.11
##    y0[13]    25.14
##    y0[14]    17.04
##    y0[15]    28.45
##    y0[16]    11.80
##    y0[17]    34.46
##    y0[18]    12.21
##    y0[19]    16.23
##    y0[20]    29.50
##    m1[1]      7.76
##    m1[2]     11.65
##    m1[3]     15.53
##    m1[4]     19.41
##    m1[5]     23.29
##    m1[6]     27.17
##    m1[7]     31.05
##    m1[8]     34.93
##    m1[9]     38.81
##    m1[10]    42.69
##    y1[1]      8.05
##    y1[2]     11.99
##    y1[3]     17.61
##    y1[4]     18.23
##    y1[5]     23.03
##    y1[6]     27.97
##    y1[7]     24.03
##    y1[8]     37.41
##    y1[9]     37.60
##    y1[10]    47.32
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -38.80 -38.37 1.49 1.06 -41.80 -37.32 1.01      370      212
##    b0       5.47   5.39 2.67 2.42   1.66   9.66 1.02      248      158
##    b1       2.06   2.07 0.27 0.24   1.63   2.47 1.02      284      307
##    s        4.75   4.60 0.87 0.76   3.61   6.40 1.00      662      726
##    m0[1]    7.68   7.62 2.41 2.20   4.23  11.47 1.02      249      161
##    m0[2]   28.46  28.48 1.21 1.15  26.50  30.35 1.00     1909     1223
##    m0[3]   29.51  29.52 1.29 1.22  27.44  31.50 1.00     1660     1196
##    m0[4]   11.19  11.18 2.00 1.86   8.29  14.39 1.02      255       72
##    m0[5]   33.95  33.98 1.69 1.58  31.16  36.57 1.00      791      862
##    m0[6]   21.12  21.09 1.13 1.08  19.32  22.94 1.00      477      488
##    m0[7]   20.53  20.50 1.16 1.09  18.72  22.40 1.01      431      484
##    m0[8]   42.79  42.83 2.69 2.43  38.33  46.87 1.01      455      464
##    m0[9]   17.24  17.20 1.39 1.27  15.15  19.46 1.02      305      280
##    m0[10]  30.78  30.81 1.39 1.30  28.55  32.90 1.00     1322      968
##    m0[11]  17.48  17.44 1.37 1.25  15.42  19.66 1.01      310      313
##    m0[12]  18.31  18.28 1.31 1.19  16.33  20.39 1.01      331      368
##    m0[13]  31.87  31.92 1.48 1.39  29.46  34.19 1.00     1066      893
##    m0[14]  17.13  17.10 1.40 1.27  15.04  19.37 1.02      303      290
##    m0[15]  30.80  30.84 1.39 1.30  28.57  32.92 1.00     1317      968
##    m0[16]  19.26  19.24 1.24 1.14  17.34  21.25 1.01      363      393
##    m0[17]  30.13  30.16 1.33 1.25  28.00  32.17 1.00     1497     1165
##    m0[18]  17.56  17.51 1.37 1.24  15.51  19.73 1.01      312      326
##    m0[19]  12.65  12.63 1.84 1.70   9.98  15.61 1.02      259       87
##    m0[20]  32.33  32.37 1.53 1.43  29.82  34.71 1.00      987      866
##    y0[1]    7.75   7.76 5.52 5.03  -1.03  16.90 1.00      865     1537
##    y0[2]   28.53  28.41 5.10 4.91  20.31  36.89 1.00     2216     2060
##    y0[3]   29.55  29.52 4.88 4.66  21.47  37.22 1.00     2120     2024
##    y0[4]   11.17  11.19 5.25 5.07   2.69  19.57 1.01      772     1239
##    y0[5]   33.91  34.01 5.13 4.71  25.59  42.23 1.00     1794     1878
##    y0[6]   21.17  21.07 4.97 4.94  13.15  29.53 1.00     1667     1624
##    y0[7]   20.63  20.52 4.84 4.68  12.80  28.78 1.00     1397     1984
##    y0[8]   42.94  42.96 5.43 5.26  34.05  51.68 1.00     1203     1372
##    y0[9]   17.26  17.20 4.91 4.87   9.43  25.31 1.00     1886     1658
##    y0[10]  30.82  30.77 5.14 5.17  22.34  39.01 1.00     1752     1816
##    y0[11]  17.46  17.55 5.05 4.92   9.15  25.73 1.00     1586     1591
##    y0[12]  18.32  18.35 4.98 4.79  10.12  26.25 1.00     1623     1805
##    y0[13]  31.74  31.80 5.09 4.85  23.29  39.61 1.00     1914     1617
##    y0[14]  17.04  17.17 5.13 4.97   8.86  25.46 1.00     1595     1728
##    y0[15]  30.83  30.97 5.04 4.95  22.58  38.66 1.00     1881     1804
##    y0[16]  19.10  19.23 5.00 4.72  10.81  27.10 1.00     1491     2055
##    y0[17]  30.19  30.33 5.02 4.90  22.07  38.22 1.00     1983     1783
##    y0[18]  17.49  17.50 5.07 4.85   9.23  25.61 1.00     1515     1665
##    y0[19]  12.75  12.66 5.30 5.17   4.16  21.31 1.00      953     1475
##    y0[20]  32.38  32.36 5.07 4.95  23.94  40.75 1.00     1868     1933
##    m1[1]    7.68   7.62 2.41 2.20   4.23  11.47 1.02      249      161
##    m1[2]   11.58  11.57 1.96 1.81   8.75  14.71 1.02      256       72
##    m1[3]   15.48  15.45 1.55 1.41  13.18  17.96 1.02      279      167
##    m1[4]   19.38  19.36 1.23 1.14  17.47  21.36 1.01      368      393
##    m1[5]   23.28  23.27 1.07 1.01  21.57  25.00 1.00      808      865
##    m1[6]   27.18  27.19 1.14 1.11  25.35  28.98 1.00     1985     1239
##    m1[7]   31.08  31.12 1.41 1.32  28.81  33.25 1.00     1241      939
##    m1[8]   34.98  35.03 1.79 1.68  32.00  37.76 1.00      707      766
##    m1[9]   38.89  38.91 2.23 2.06  35.13  42.31 1.01      532      493
##    m1[10]  42.79  42.83 2.69 2.43  38.33  46.87 1.01      455      464
##    y1[1]    7.73   7.82 5.38 5.12  -0.99  16.36 1.00      921      937
##    y1[2]   11.47  11.39 5.24 4.98   2.87  20.28 1.00      923     1216
##    y1[3]   15.24  15.32 5.14 4.92   6.60  23.72 1.00     1058     1339
##    y1[4]   19.76  19.85 5.01 4.74  11.46  27.88 1.00     1502     1687
##    y1[5]   23.15  23.35 4.86 4.81  15.26  30.80 1.00     1919     1931
##    y1[6]   27.24  27.27 4.83 4.53  19.58  35.03 1.00     1811     1842
##    y1[7]   31.13  30.98 5.02 4.59  22.65  39.90 1.00     1893     1883
##    y1[8]   34.78  34.88 5.17 4.93  26.14  43.30 1.00     1972     1760
##    y1[9]   38.88  38.95 5.21 4.94  30.21  47.34 1.01     1129     1359
##    y1[10]  42.70  42.74 5.56 5.27  33.46  51.56 1.00      947     1115

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.4 seconds.
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 4 finished in 0.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -44.95 -46.14 10.26 8.76 -59.90 -25.84 1.04      118      181
##    b0       5.83   5.92  2.58 2.50   1.69  10.02 1.00      749      993
##    b1       2.02   2.02  0.27 0.26   1.58   2.46 1.00      675      738
##    s        3.58   3.63  1.37 1.38   1.36   5.77 1.03      175      130
##    sx       1.46   1.43  0.82 0.89   0.31   2.84 1.05       70      171
##    x0[1]    1.02   1.07  1.17 0.92  -1.09   2.91 1.02     1280      803
##    x0[2]   12.02  11.83  1.28 1.15  10.23  14.33 1.01      547      922
##    x0[3]   13.51  13.17  1.85 1.96  11.20  16.93 1.02      265      639
##    x0[4]    3.55   3.42  1.20 1.13   1.78   5.60 1.01      608     2135
##    x0[5]   13.38  13.46  1.14 0.97  11.53  15.12 1.01     1097     1684
##    x0[6]    6.96   7.11  1.14 1.01   4.96   8.58 1.01      693     1061
##    x0[7]    6.30   6.45  1.31 1.26   4.05   8.13 1.02      456      913
##    x0[8]   19.05  18.70  1.52 1.22  17.17  22.08 1.01      524      535
##    x0[9]    4.92   5.11  1.25 1.10   2.61   6.68 1.01      516      861
##    x0[10]  12.10  12.13  1.05 0.86  10.41  13.76 1.01     1459     1482
##    x0[11]   5.76   5.81  1.06 0.81   3.96   7.38 1.01     2777     1361
##    x0[12]   6.80   6.67  1.12 1.01   5.14   8.77 1.01      851     1908
##    x0[13]  13.88  13.62  1.39 1.27  12.08  16.52 1.02      421      870
##    x0[14]   5.92   5.88  1.05 0.90   4.25   7.64 1.01     1909     1918
##    x0[15]  12.06  12.10  1.09 0.89  10.20  13.90 1.01     2147     1661
##    x0[16]   6.42   6.52  1.04 0.84   4.58   8.03 1.01     2140     1140
##    x0[17]  10.39  10.60  1.56 1.72   7.70  12.44 1.03      177     1069
##    x0[18]   6.80   6.64  1.24 1.23   5.04   8.89 1.01      482     1396
##    x0[19]   2.82   3.03  1.25 0.99   0.54   4.57 1.01      736      897
##    x0[20]  11.59  11.83  1.47 1.60   9.02  13.54 1.02      268     1117
##    m0[1]    8.01   8.03  2.48 2.30   3.89  12.04 1.00     2213     2134
##    m0[2]   30.06  30.02  2.46 2.52  26.17  34.13 1.01      613     1630
##    m0[3]   33.02  32.74  3.46 4.05  27.88  38.66 1.02      269     1754
##    m0[4]   13.08  13.11  2.59 2.64   8.84  17.19 1.01      517     1923
##    m0[5]   32.82  32.76  2.40 2.32  29.02  36.89 1.01     1271     2577
##    m0[6]   19.95  20.03  2.26 2.22  16.23  23.46 1.01      814     2204
##    m0[7]   18.63  18.74  2.57 2.70  14.37  22.62 1.01      474     1729
##    m0[8]   44.16  44.35  2.91 2.86  39.12  48.63 1.01      760     2547
##    m0[9]   15.85  15.91  2.47 2.52  11.90  19.80 1.01      666     3058
##    m0[10]  30.25  30.21  2.15 1.93  26.79  33.94 1.00     1964     2530
##    m0[11]  17.53  17.56  2.07 1.84  14.04  20.95 1.00     3322     2484
##    m0[12]  19.59  19.58  2.29 2.29  16.00  23.25 1.01      747     2154
##    m0[13]  33.78  33.77  2.63 2.82  29.59  38.02 1.01      436     2279
##    m0[14]  17.83  17.88  2.15 2.02  14.24  21.31 1.00     1835     2376
##    m0[15]  30.18  30.17  2.26 2.11  26.52  33.88 1.00     3026     2206
##    m0[16]  18.86  18.87  2.09 1.87  15.37  22.18 1.00     3041     2615
##    m0[17]  26.85  27.07  3.22 3.64  21.64  31.71 1.03      165     1437
##    m0[18]  19.58  19.50  2.54 2.68  15.54  23.67 1.02      463     2091
##    m0[19]  11.63  11.57  2.44 2.40   7.65  15.68 1.00     1064     2328
##    m0[20]  29.27  29.37  3.13 3.52  24.23  34.14 1.02      279     1571
##    y0[1]    8.06   8.02  4.54 4.00   0.72  15.57 1.00     2867     2468
##    y0[2]   30.03  30.43  4.59 4.19  21.95  36.85 1.00     1947     3122
##    y0[3]   33.02  33.42  5.19 5.18  24.02  40.59 1.02      444     2752
##    y0[4]   13.04  13.38  4.63 4.12   4.90  20.25 1.00     1877     2511
##    y0[5]   32.93  32.62  4.44 3.87  25.96  40.65 1.00     2237     3209
##    y0[6]   19.96  19.67  4.51 4.06  12.75  27.88 1.00     2369     2697
##    y0[7]   18.62  18.14  4.66 4.32  11.78  26.77 1.00     1626     3471
##    y0[8]   44.22  44.57  4.79 4.29  36.06  51.48 1.00     2031     2766
##    y0[9]   15.87  15.50  4.54 4.10   9.00  23.65 1.00     1788     3538
##    y0[10]  30.25  30.24  4.39 3.90  22.99  37.41 1.00     3540     3092
##    y0[11]  17.57  17.55  4.34 3.86  10.31  24.83 1.00     3680     3348
##    y0[12]  19.55  19.80  4.49 3.95  11.80  26.47 1.00     2457     3358
##    y0[13]  33.90  34.31  4.54 4.23  25.78  40.50 1.00     1507     2918
##    y0[14]  17.85  18.02  4.37 3.83  10.47  24.71 1.00     3754     3301
##    y0[15]  30.13  29.98  4.34 3.83  23.12  37.43 1.00     4201     2938
##    y0[16]  18.75  18.65  4.36 3.85  11.64  25.97 1.00     4073     3483
##    y0[17]  26.75  26.22  4.95 4.83  19.62  35.48 1.01      692     1964
##    y0[18]  19.58  20.03  4.61 4.37  11.63  26.47 1.00     1569     2643
##    y0[19]  11.61  11.30  4.54 4.14   4.55  19.26 1.00     2653     3083
##    y0[20]  29.23  28.78  4.96 4.86  22.14  38.09 1.01      665     2501
##    m1[1]    8.00   8.05  2.33 2.25   4.30  11.77 1.00      771     1011
##    m1[2]   11.82  11.84  1.90 1.81   8.77  14.93 1.00      838     1231
##    m1[3]   15.64  15.67  1.53 1.45  13.17  18.15 1.00      994     1394
##    m1[4]   19.47  19.46  1.25 1.17  17.44  21.49 1.00     1354     1500
##    m1[5]   23.29  23.30  1.13 1.08  21.44  25.14 1.00     1797     1822
##    m1[6]   27.11  27.12  1.24 1.19  25.08  29.15 1.00     1432     2294
##    m1[7]   30.94  30.96  1.51 1.44  28.44  33.41 1.00     1060     1840
##    m1[8]   34.76  34.76  1.88 1.81  31.66  37.84 1.00      891     1468
##    m1[9]   38.58  38.56  2.31 2.23  34.80  42.34 1.00      806     1186
##    m1[10]  42.41  42.39  2.76 2.68  37.88  46.87 1.00      763     1018
##    x1[1]    1.07   1.09  1.67 1.18  -1.66   3.73 1.02     3843     2173
##    x1[2]    2.97   2.96  1.69 1.14   0.20   5.81 1.02     3750     1971
##    x1[3]    4.87   4.85  1.68 1.15   2.17   7.64 1.01     3855     2675
##    x1[4]    6.73   6.73  1.70 1.18   4.06   9.41 1.02     3709     2192
##    x1[5]    8.63   8.62  1.65 1.14   5.90  11.28 1.02     4148     2087
##    x1[6]   10.53  10.53  1.66 1.14   7.77  13.28 1.02     4013     2495
##    x1[7]   12.46  12.45  1.68 1.18   9.80  15.16 1.02     4024     1813
##    x1[8]   14.35  14.34  1.66 1.13  11.59  16.95 1.02     4028     2139
##    x1[9]   16.23  16.20  1.68 1.18  13.57  19.08 1.02     3892     2389
##    x1[10]  18.11  18.12  1.67 1.14  15.46  20.83 1.02     3962     2205
##    y1[1]    7.95   7.80  4.59 4.13   0.62  15.38 1.00     1619     2319
##    y1[2]   11.82  11.78  4.34 3.88   4.60  18.92 1.00     2270     1950
##    y1[3]   15.68  15.71  4.16 3.60   9.03  22.43 1.00     3144     3702
##    y1[4]   19.45  19.45  3.88 3.40  12.96  25.86 1.00     3779     3648
##    y1[5]   23.20  23.28  3.99 3.43  16.55  29.73 1.00     3594     2668
##    y1[6]   26.99  27.03  4.11 3.63  20.31  33.64 1.00     3395     3368
##    y1[7]   30.93  30.91  4.10 3.61  24.11  37.57 1.00     3166     2843
##    y1[8]   34.76  34.80  4.23 3.78  27.89  41.60 1.00     2548     2759
##    y1[9]   38.75  38.71  4.49 4.10  31.52  45.92 1.00     1714     2471
##    y1[10]  42.44  42.44  4.60 4.26  34.90  50.07 1.00     1454     2259

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -8140.49 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20       -32.945   0.000939809   2.02029e-05           1           1       39    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -32.94
##    b0         5.88
##    b1         1.97
##    s          3.15
##    m0[1]     11.78
##    m0[2]     10.27
##    m0[3]     20.62
##    m0[4]     21.91
##    m0[5]     13.05
##    m0[6]      6.70
##    m0[7]     10.94
##    m0[8]     18.87
##    m0[9]     12.98
##    m0[10]    20.66
##    m0[11]    22.87
##    m0[12]     6.77
##    m0[13]    13.26
##    m0[14]    16.58
##    m0[15]    23.05
##    m0[16]    13.73
##    m0[17]    12.55
##    m0[18]     6.07
##    m0[19]    17.99
##    m0[20]    11.58
##    y0[1]      9.84
##    y0[2]      9.29
##    y0[3]     18.63
##    y0[4]     18.17
##    y0[5]     13.89
##    y0[6]      8.13
##    y0[7]     16.62
##    y0[8]     14.74
##    y0[9]     11.95
##    y0[10]    23.93
##    y0[11]    26.61
##    y0[12]     4.81
##    y0[13]    10.59
##    y0[14]    16.98
##    y0[15]    32.01
##    y0[16]    12.84
##    y0[17]    11.85
##    y0[18]     7.09
##    y0[19]     8.88
##    y0[20]    10.60
##    m1[1]      6.07
##    m1[2]      7.95
##    m1[3]      9.84
##    m1[4]     11.73
##    m1[5]     13.61
##    m1[6]     15.50
##    m1[7]     17.39
##    m1[8]     19.27
##    m1[9]     21.16
##    m1[10]    23.05
##    y1[1]     12.71
##    y1[2]      5.00
##    y1[3]      6.60
##    y1[4]     15.00
##    y1[5]     13.39
##    y1[6]     15.26
##    y1[7]     17.53
##    y1[8]     20.62
##    y1[9]     24.12
##    y1[10]    25.35
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.47 -33.16 1.34 1.16 -35.96 -31.96 1.01      529      714
##    b0       5.81   5.80 1.67 1.57   2.95   8.57 1.01      299      283
##    b1       1.98   1.99 0.31 0.30   1.46   2.51 1.01      393      604
##    s        3.57   3.49 0.65 0.60   2.73   4.72 1.00     1470     1506
##    m0[1]   11.75  11.72 0.98 0.96  10.17  13.40 1.01      387      477
##    m0[2]   10.23  10.22 1.12 1.07   8.43  12.12 1.01      334      322
##    m0[3]   20.65  20.66 1.24 1.19  18.61  22.68 1.00     1549     1383
##    m0[4]   21.95  21.96 1.40 1.34  19.67  24.25 1.00     1225     1308
##    m0[5]   13.03  13.01 0.89 0.87  11.58  14.52 1.00      484      604
##    m0[6]    6.63   6.63 1.56 1.48   4.01   9.24 1.01      300      276
##    m0[7]   10.90  10.88 1.06 1.01   9.19  12.69 1.01      351      360
##    m0[8]   18.89  18.88 1.05 1.00  17.12  20.58 1.00     1905     1448
##    m0[9]   12.95  12.93 0.90 0.88  11.50  14.45 1.00      475      600
##    m0[10]  20.69  20.70 1.24 1.19  18.64  22.72 1.00     1541     1383
##    m0[11]  22.91  22.94 1.52 1.45  20.43  25.40 1.00     1027     1290
##    m0[12]   6.70   6.71 1.55 1.47   4.10   9.29 1.01      300      282
##    m0[13]  13.24  13.22 0.88 0.86  11.82  14.71 1.00      511      635
##    m0[14]  16.59  16.59 0.88 0.88  15.10  18.02 1.00     1760     1513
##    m0[15]  23.09  23.12 1.55 1.47  20.57  25.61 1.00      995     1290
##    m0[16]  13.72  13.70 0.86 0.86  12.34  15.14 1.00      580      649
##    m0[17]  12.53  12.50 0.92 0.90  11.02  14.08 1.00      436      534
##    m0[18]   5.99   5.99 1.65 1.55   3.19   8.73 1.01      299      274
##    m0[19]  18.01  17.99 0.97 0.95  16.38  19.58 1.00     1957     1542
##    m0[20]  11.55  11.52 1.00 0.97   9.94  13.23 1.01      377      461
##    y0[1]   11.77  11.77 3.80 3.62   5.46  18.03 1.00     1489     1826
##    y0[2]   10.20  10.19 3.81 3.62   4.07  16.37 1.00     1256     1692
##    y0[3]   20.53  20.60 3.83 3.77  14.22  26.55 1.00     1954     2012
##    y0[4]   21.98  22.05 3.91 3.85  15.63  28.13 1.00     1746     1972
##    y0[5]   13.01  13.03 3.72 3.63   7.06  19.00 1.00     1752     1794
##    y0[6]    6.59   6.62 4.11 3.96  -0.48  13.03 1.00      947     1720
##    y0[7]   10.99  11.01 3.75 3.62   4.84  17.18 1.00     1934     1382
##    y0[8]   18.87  18.96 3.76 3.70  12.59  25.24 1.00     2068     1905
##    y0[9]   12.94  12.99 3.73 3.58   6.61  18.96 1.00     1930     1997
##    y0[10]  20.71  20.69 3.94 3.66  14.12  27.19 1.00     1845     1675
##    y0[11]  22.85  22.83 4.08 3.75  16.36  29.65 1.00     1691     1483
##    y0[12]   6.67   6.82 3.94 3.81   0.18  13.14 1.00     1201     1475
##    y0[13]  13.24  13.29 3.67 3.56   7.23  19.13 1.00     1914     1968
##    y0[14]  16.62  16.57 3.72 3.67  10.53  22.54 1.00     1799     1893
##    y0[15]  23.04  22.99 4.01 3.92  16.53  29.69 1.00     1943     1893
##    y0[16]  13.54  13.54 3.71 3.53   7.58  19.58 1.00     1914     1732
##    y0[17]  12.49  12.50 3.75 3.68   6.37  18.78 1.00     1713     1587
##    y0[18]   5.89   6.02 4.08 3.76  -1.03  12.10 1.00     1038     1407
##    y0[19]  18.05  18.03 3.77 3.58  11.97  24.11 1.00     2095     1972
##    y0[20]  11.47  11.54 3.81 3.66   5.23  17.59 1.00     1556     1650
##    m1[1]    5.99   5.99 1.65 1.55   3.19   8.73 1.01      299      274
##    m1[2]    7.89   7.91 1.40 1.33   5.60  10.25 1.01      306      275
##    m1[3]    9.79   9.79 1.17 1.12   7.94  11.74 1.01      326      314
##    m1[4]   11.69  11.66 0.98 0.97  10.10  13.36 1.01      384      477
##    m1[5]   13.59  13.58 0.87 0.86  12.20  15.03 1.00      561      634
##    m1[6]   15.49  15.50 0.85 0.83  14.14  16.87 1.00     1091     1140
##    m1[7]   17.39  17.39 0.93 0.92  15.81  18.90 1.00     1916     1554
##    m1[8]   19.29  19.28 1.09 1.03  17.45  21.06 1.00     1869     1518
##    m1[9]   21.19  21.20 1.30 1.25  19.06  23.34 1.00     1440     1364
##    m1[10]  23.09  23.12 1.55 1.47  20.57  25.61 1.00      995     1290
##    y1[1]    5.87   5.81 4.00 3.88  -0.54  12.36 1.00     1204     1512
##    y1[2]    7.95   7.94 3.86 3.85   1.88  14.40 1.00      970     1846
##    y1[3]    9.78   9.82 3.76 3.63   3.48  15.98 1.00     1561     1755
##    y1[4]   11.63  11.70 3.83 3.80   5.34  17.81 1.00     1763     1658
##    y1[5]   13.70  13.67 3.78 3.59   7.69  19.94 1.00     1812     1892
##    y1[6]   15.39  15.35 3.67 3.48   9.52  21.54 1.00     1891     1931
##    y1[7]   17.27  17.30 3.64 3.52  11.35  23.24 1.00     1797     1880
##    y1[8]   19.30  19.20 3.65 3.49  13.46  25.33 1.00     2026     2013
##    y1[9]   21.25  21.32 3.92 3.81  14.84  27.58 1.00     1943     1822
##    y1[10]  23.20  23.23 3.98 3.70  16.82  29.83 1.00     2004     1839

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -124.542 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24      -9.48434   0.000122605    0.00219896      0.9709      0.9709       34    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -9.48
##    b0         4.45
##    b1         2.07
##    s          0.47
##    m0[1]     10.66
##    m0[2]      9.07
##    m0[3]     19.95
##    m0[4]     21.30
##    m0[5]     11.99
##    m0[6]      5.31
##    m0[7]      9.77
##    m0[8]     18.11
##    m0[9]     11.92
##    m0[10]    19.99
##    m0[11]    22.31
##    m0[12]     5.39
##    m0[13]    12.22
##    m0[14]    15.71
##    m0[15]    22.50
##    m0[16]    12.71
##    m0[17]    11.47
##    m0[18]     4.65
##    m0[19]    17.19
##    m0[20]    10.45
##    y0[1]     11.60
##    y0[2]      8.74
##    y0[3]     21.60
##    y0[4]     21.65
##    y0[5]     11.81
##    y0[6]      8.40
##    y0[7]      9.61
##    y0[8]     18.20
##    y0[9]     12.06
##    y0[10]    19.18
##    y0[11]    19.53
##    y0[12]     4.47
##    y0[13]    12.48
##    y0[14]    15.24
##    y0[15]    21.71
##    y0[16]    12.76
##    y0[17]    11.65
##    y0[18]     4.39
##    y0[19]    19.67
##    y0[20]     9.57
##    m1[1]      4.65
##    m1[2]      6.63
##    m1[3]      8.62
##    m1[4]     10.60
##    m1[5]     12.58
##    m1[6]     14.57
##    m1[7]     16.55
##    m1[8]     18.54
##    m1[9]     20.52
##    m1[10]    22.50
##    y1[1]      4.61
##    y1[2]      7.05
##    y1[3]      8.36
##    y1[4]     12.84
##    y1[5]     12.42
##    y1[6]     15.57
##    y1[7]     17.02
##    y1[8]     25.72
##    y1[9]     19.62
##    y1[10]    20.69
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -11.99 -11.61   1.44 1.16 -15.02 -10.42 1.01      531      594
##    b0       4.45   4.45   0.37 0.33   3.84   5.07 1.00      669      626
##    b1       2.08   2.08   0.07 0.06   1.98   2.20 1.00      777     1001
##    s        0.60   0.58   0.20 0.19   0.34   0.97 1.01      525      286
##    m0[1]   10.70  10.69   0.24 0.23  10.33  11.11 1.01      853      766
##    m0[2]    9.11   9.10   0.26 0.25   8.69   9.55 1.00      759      685
##    m0[3]   20.07  20.03   0.30 0.27  19.64  20.63 1.00     1702     1380
##    m0[4]   21.43  21.39   0.33 0.29  20.96  22.05 1.00     1540     1299
##    m0[5]   12.05  12.03   0.22 0.22  11.71  12.45 1.01      995     1003
##    m0[6]    5.32   5.31   0.35 0.32   4.74   5.90 1.00      675      621
##    m0[7]    9.81   9.80   0.25 0.24   9.42  10.24 1.00      794      699
##    m0[8]   18.21  18.19   0.26 0.23  17.84  18.70 1.00     1764     1378
##    m0[9]   11.97  11.96   0.22 0.22  11.63  12.37 1.01      984     1003
##    m0[10]  20.11  20.07   0.30 0.27  19.68  20.67 1.00     1696     1380
##    m0[11]  22.45  22.41   0.36 0.32  21.94  23.11 1.00     1440     1253
##    m0[12]   5.40   5.39   0.34 0.32   4.82   5.97 1.00      675      621
##    m0[13]  12.28  12.26   0.22 0.21  11.94  12.67 1.01     1026     1047
##    m0[14]  15.79  15.77   0.23 0.21  15.47  16.20 1.00     1656     1287
##    m0[15]  22.64  22.60   0.36 0.32  22.12  23.31 1.00     1422     1258
##    m0[16]  12.77  12.75   0.22 0.21  12.44  13.16 1.01     1103     1207
##    m0[17]  11.52  11.51   0.23 0.22  11.17  11.93 1.01      932      896
##    m0[18]   4.65   4.65   0.36 0.33   4.05   5.26 1.00      670      626
##    m0[19]  17.29  17.26   0.25 0.22  16.94  17.74 1.00     1755     1436
##    m0[20]  10.49  10.48   0.24 0.23  10.11  10.91 1.01      837      709
##    y0[1]   10.60  10.68  20.96 0.90   6.13  14.42 1.00     1809     2015
##    y0[2]    9.37   9.10   9.44 0.98   5.27  13.39 1.00     1915     1838
##    y0[3]   19.68  20.05  31.14 1.01  15.90  24.55 1.00     1971     1892
##    y0[4]   26.33  21.44 193.08 0.93  18.01  25.25 1.00     1930     1793
##    y0[5]   12.50  12.03  19.59 0.91   8.72  16.22 1.00     2048     1931
##    y0[6]    5.46   5.31  10.35 1.03   1.64   9.50 1.00     1704     1787
##    y0[7]    9.28   9.82  31.06 0.95   5.65  13.58 1.00     1936     1843
##    y0[8]   18.69  18.16  16.42 0.90  14.64  22.03 1.00     2067     1877
##    y0[9]   12.28  11.94  20.98 0.94   8.13  15.61 1.00     1795     1973
##    y0[10]  19.98  20.07  17.56 0.97  15.56  23.94 1.00     1901     1974
##    y0[11]  22.25  22.40   8.54 0.98  18.75  26.16 1.00     1931     1857
##    y0[12]   2.80   5.42  79.01 0.98   1.74   9.30 1.00     1868     2052
##    y0[13]  12.18  12.27   9.12 0.96   8.38  15.68 1.00     1980     1980
##    y0[14]  15.82  15.77  10.70 0.85  11.99  19.71 1.00     2129     2036
##    y0[15]  23.04  22.60  18.71 0.95  18.65  26.72 1.00     1975     1822
##    y0[16]  12.75  12.79  45.54 0.87   8.75  16.04 1.00     1995     1694
##    y0[17]  11.91  11.54  23.89 0.92   7.96  15.45 1.00     2089     1835
##    y0[18]   4.36   4.65  23.74 1.00   0.73   8.31 1.00     1622     1755
##    y0[19]  16.60  17.25  36.10 0.85  13.46  21.05 1.00     2140     1972
##    y0[20]  10.63  10.49   7.08 0.90   7.49  14.11 1.00     1927     1787
##    m1[1]    4.65   4.65   0.36 0.33   4.05   5.26 1.00      670      626
##    m1[2]    6.65   6.64   0.31 0.29   6.13   7.17 1.00      691      612
##    m1[3]    8.65   8.64   0.27 0.26   8.22   9.10 1.00      740      685
##    m1[4]   10.65  10.63   0.24 0.23  10.27  11.06 1.01      849      766
##    m1[5]   12.65  12.63   0.22 0.21  12.31  13.03 1.01     1081     1137
##    m1[6]   14.64  14.62   0.22 0.21  14.32  15.03 1.00     1458     1342
##    m1[7]   16.64  16.61   0.24 0.22  16.31  17.08 1.00     1721     1323
##    m1[8]   18.64  18.61   0.27 0.24  18.26  19.14 1.00     1762     1319
##    m1[9]   20.64  20.60   0.31 0.28  20.20  21.22 1.00     1630     1322
##    m1[10]  22.64  22.60   0.36 0.32  22.12  23.31 1.00     1422     1258
##    y1[1]    5.24   4.61  28.58 0.94   0.41   7.86 1.00     1578     1631
##    y1[2]    6.35   6.64  19.10 0.91   2.88  10.33 1.00     1602     1812
##    y1[3]    9.78   8.64  24.87 0.97   5.15  12.72 1.00     1801     1987
##    y1[4]   10.90  10.66  34.15 0.94   6.65  14.53 1.00     1942     1933
##    y1[5]   12.92  12.63  19.90 0.90   8.44  16.48 1.00     1945     1894
##    y1[6]   15.50  14.65  49.43 0.93  10.52  18.58 1.00     1947     1896
##    y1[7]   16.40  16.64  16.06 0.88  12.98  20.80 1.00     2082     1707
##    y1[8]   17.40  18.62  29.80 0.91  14.62  22.55 1.00     1867     1784
##    y1[9]   20.66  20.62  12.96 0.95  16.71  24.68 1.00     1924     1969
##    y1[10]  17.45  22.62 193.75 1.00  19.09  26.13 1.00     2097     1755

censored data

objective variable has NA

(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)

ex11-0.stan

data{
  int N0;
  array[N0] vector[2] xy;
  int N1;
  vector[N1] x1;
}
parameters{
  vector[2] m;
  cov_matrix[2] s;
}
model{
  target+=multi_normal_lpdf(xy | m, s);
  x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
  vector[2] xy1;
  xy1=multi_normal_rng(m,s);
  real cr;
  cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}
n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)
## [1] 0.905
qplot(x,y)

L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)
## [1] 0.716
qplot(x0,y0)

mdl=cmdstan_model('./ex11-0.stan') 

data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -107248 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       59      -106.626   0.000515973     0.0018835           1           1       79    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    -106.63
##    m[1]       4.61
##    m[2]      25.09
##    s[1,1]     7.31
##    s[2,1]    20.79
##    s[1,2]    20.79
##    s[2,2]    72.51
##    xy1[1]     5.82
##    xy1[2]    25.84
##    cr         0.90
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.5 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -102.16 -101.85  1.62  1.46 -105.22 -100.13 1.00      664     1041
##    m[1]      4.61    4.62  0.54  0.50    3.68    5.49 1.00      726      708
##    m[2]     25.16   25.26  2.56  2.35   20.90   29.31 1.01      400      374
##    s[1,1]    9.43    8.92  2.83  2.50    5.80   14.79 1.01     1034     1347
##    s[2,1]   26.35   24.53 10.92  9.61   11.86   46.80 1.01      437      608
##    s[1,2]   26.35   24.53 10.92  9.61   11.86   46.80 1.01      437      608
##    s[2,2]  100.27   88.62 49.64 40.97   41.46  198.19 1.01      425      517
##    xy1[1]    4.60    4.61  3.14  2.99   -0.61    9.70 1.00     1631     1673
##    xy1[2]   25.27   25.60 10.35  9.43    8.40   41.68 1.00     1468     1346
##    cr        0.86    0.89  0.12  0.07    0.66    0.96 1.01      529      630

xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])
## [1] 0.835
qplot(xy[,,1],xy[,,2])

objective variable is censored

y i=1-N, y~N(m,s)  
  actual          ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
mdl=cmdstan_model('./ex8-0.stan')

data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -19025.6 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
##       24      -7.43161   0.000242549    0.00045938      0.9732      0.9732       53    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -7.43
##    b0        15.10
##    b1         1.94
##    s          1.39
##    m0[1]     20.06
##    m0[2]     26.77
##    m0[3]     28.06
##    m0[4]     24.88
##    m0[5]     18.24
##    m0[6]     19.16
##    m0[7]     28.15
##    m0[8]     20.52
##    m0[9]     19.25
##    y0[1]     20.67
##    y0[2]     26.52
##    y0[3]     30.90
##    y0[4]     26.40
##    y0[5]     19.86
##    y0[6]     20.52
##    y0[7]     27.71
##    y0[8]     20.59
##    y0[9]     18.69
##    m1[1]     15.79
##    m1[2]     17.61
##    m1[3]     19.42
##    m1[4]     21.24
##    m1[5]     23.05
##    m1[6]     24.87
##    m1[7]     26.69
##    m1[8]     28.50
##    m1[9]     30.32
##    m1[10]    32.13
##    y1[1]     15.53
##    y1[2]     17.54
##    y1[3]     20.63
##    y1[4]     21.65
##    y1[5]     24.75
##    y1[6]     26.84
##    y1[7]     26.62
##    y1[8]     29.61
##    y1[9]     29.45
##    y1[10]    32.79
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -8.84  -8.48 1.48 1.19 -11.78 -7.27 1.00      446      310
##    b0     15.17  15.16 1.48 1.32  12.89 17.49 1.02      270      315
##    b1      1.93   1.93 0.33 0.29   1.40  2.46 1.02      342      463
##    s       1.93   1.79 0.68 0.53   1.19  3.14 1.01      605      526
##    m0[1]  20.11  20.09 0.83 0.70  18.86 21.47 1.02      335      311
##    m0[2]  26.78  26.77 0.98 0.87  25.30 28.35 1.00     1491      789
##    m0[3]  28.06  28.05 1.14 1.01  26.30 29.89 1.00     1126      684
##    m0[4]  24.90  24.90 0.79 0.67  23.72 26.16 1.00     1859      759
##    m0[5]  18.29  18.27 1.04 0.90  16.70 19.91 1.02      279      327
##    m0[6]  19.21  19.19 0.93 0.80  17.81 20.65 1.02      292      333
##    m0[7]  28.16  28.14 1.16 1.02  26.37 30.00 1.00     1106      685
##    m0[8]  20.57  20.54 0.79 0.66  19.37 21.86 1.01      382      388
##    m0[9]  19.30  19.27 0.92 0.79  17.91 20.72 1.02      293      335
##    y0[1]  20.19  20.15 2.19 1.83  16.73 23.70 1.00     1677     1222
##    y0[2]  26.79  26.76 2.31 1.96  23.11 30.36 1.00     1562     1234
##    y0[3]  28.04  28.06 2.28 2.05  24.30 31.68 1.00     1718     1417
##    y0[4]  24.87  24.81 2.21 1.90  21.45 28.32 1.00     1955     1226
##    y0[5]  18.21  18.24 2.29 1.99  14.61 21.75 1.00     1145     1296
##    y0[6]  19.30  19.26 2.32 1.91  15.71 22.75 1.00      960     1186
##    y0[7]  28.11  28.15 2.36 2.05  24.08 31.79 1.00     1761     1595
##    y0[8]  20.59  20.54 2.26 1.79  17.26 23.97 1.00     1279     1188
##    y0[9]  19.20  19.20 2.19 1.94  15.63 22.61 1.00     1122     1257
##    m1[1]  15.86  15.84 1.38 1.21  13.76 18.02 1.02      270      316
##    m1[2]  17.67  17.64 1.12 0.97  15.93 19.40 1.02      274      321
##    m1[3]  19.47  19.45 0.90 0.78  18.11 20.88 1.02      297      335
##    m1[4]  21.28  21.22 0.74 0.60  20.18 22.51 1.01      461      382
##    m1[5]  23.08  23.06 0.70 0.58  22.06 24.20 1.00     1072      552
##    m1[6]  24.89  24.88 0.78 0.67  23.71 26.15 1.00     1860      759
##    m1[7]  26.70  26.69 0.97 0.86  25.23 28.24 1.01     1510      795
##    m1[8]  28.50  28.49 1.20 1.06  26.63 30.42 1.01     1028      703
##    m1[9]  30.31  30.29 1.47 1.26  27.99 32.65 1.01      788      580
##    m1[10] 32.11  32.09 1.75 1.51  29.39 34.90 1.01      668      636
##    y1[1]  15.77  15.71 2.48 2.14  11.80 19.51 1.01      718      934
##    y1[2]  17.63  17.62 2.34 2.01  13.98 21.31 1.01      990     1194
##    y1[3]  19.54  19.46 2.28 1.97  15.97 23.22 1.00     1292     1602
##    y1[4]  21.33  21.32 2.19 1.92  17.98 24.74 1.00     1697     1700
##    y1[5]  23.10  23.10 2.14 1.89  19.73 26.63 1.00     1839     1656
##    y1[6]  24.83  24.82 2.27 1.89  21.45 28.08 1.00     1758     1817
##    y1[7]  26.61  26.58 2.32 1.96  22.98 30.21 1.00     1839     1662
##    y1[8]  28.52  28.49 2.46 2.06  24.69 32.42 1.00     1790     1277
##    y1[9]  30.30  30.35 2.42 2.09  26.47 34.01 1.00     1889     1485
##    y1[10] 32.17  32.20 2.84 2.40  27.89 36.41 1.00     1161      946

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan') 

mle=mdl$optimize(data=data)
## Initial log joint probability = -205.788 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       31      -16.1894   8.93464e-05     0.0002329        0.88        0.88       42    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -16.19
##    b0        11.35
##    b1         2.74
##    s          2.27
##    m0[1]     18.37
##    m0[2]     27.85
##    m0[3]     29.68
##    m0[4]     25.19
##    m0[5]     15.79
##    m0[6]     17.09
##    m0[7]     29.81
##    m0[8]     19.02
##    m0[9]     17.22
##    y0[1]     17.99
##    y0[2]     27.24
##    y0[3]     28.81
##    y0[4]     26.97
##    y0[5]     14.17
##    y0[6]     19.78
##    y0[7]     29.48
##    y0[8]     20.25
##    y0[9]     19.30
##    m1[1]     12.33
##    m1[2]     14.89
##    m1[3]     17.46
##    m1[4]     20.03
##    m1[5]     22.60
##    m1[6]     25.17
##    m1[7]     27.73
##    m1[8]     30.30
##    m1[9]     32.87
##    m1[10]    35.44
##    y1[1]      9.00
##    y1[2]     10.57
##    y1[3]     21.79
##    y1[4]     21.48
##    y1[5]     22.23
##    y1[6]     24.70
##    y1[7]     28.68
##    y1[8]     31.69
##    y1[9]     32.46
##    y1[10]    37.90
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -17.19 -16.74 1.71 1.15 -20.16 -15.53 1.02      218      156
##    b0       9.94  10.43 2.89 1.70   5.98  12.92 1.03      132      113
##    b1       3.10   2.98 0.71 0.40   2.42   4.05 1.02      159      163
##    s        3.31   2.97 1.66 0.85   1.90   5.60 1.02      229      208
##    m0[1]   17.87  18.08 1.37 0.97  15.69  19.52 1.03      173      151
##    m0[2]   28.59  28.33 1.92 1.21  26.60  31.20 1.01      810      233
##    m0[3]   30.65  30.30 2.33 1.43  28.38  33.69 1.01      482      216
##    m0[4]   25.58  25.45 1.40 0.93  23.93  27.45 1.00     1409      489
##    m0[5]   14.96  15.27 1.87 1.21  12.11  16.97 1.03      140      138
##    m0[6]   16.43  16.70 1.60 1.09  13.92  18.23 1.03      145      143
##    m0[7]   30.80  30.44 2.36 1.45  28.51  33.89 1.01      469      216
##    m0[8]   18.61  18.77 1.27 0.93  16.65  20.18 1.03      207      201
##    m0[9]   16.57  16.84 1.58 1.08  14.10  18.34 1.03      145      144
##    y0[1]   17.84  18.10 3.96 3.11  11.90  23.21 1.00     1152      595
##    y0[2]   28.63  28.38 4.13 3.29  23.00  35.11 1.01     1714      814
##    y0[3]   30.59  30.20 4.50 3.43  24.60  37.29 1.01     1581      566
##    y0[4]   25.56  25.41 3.78 3.25  19.67  31.65 1.00     1830      852
##    y0[5]   14.93  15.24 4.33 3.21   8.51  20.65 1.01      754      278
##    y0[6]   16.45  16.63 3.92 3.14  10.30  21.76 1.01     1117      279
##    y0[7]   30.77  30.51 4.24 3.33  24.73  37.45 1.00     1620      538
##    y0[8]   18.69  18.70 4.03 3.07  12.60  24.43 1.00     1787      780
##    y0[9]   16.65  16.87 4.19 3.08  10.45  22.55 1.01      905      377
##    m1[1]   11.04  11.50 2.66 1.60   7.39  13.81 1.03      132      122
##    m1[2]   13.95  14.29 2.06 1.32  10.90  16.14 1.03      137      140
##    m1[3]   16.85  17.10 1.53 1.05  14.46  18.59 1.03      146      146
##    m1[4]   19.75  19.87 1.15 0.86  17.94  21.26 1.02      363      210
##    m1[5]   22.65  22.65 1.09 0.82  21.10  24.17 1.01     1103      712
##    m1[6]   25.55  25.43 1.39 0.93  23.91  27.43 1.00     1415      500
##    m1[7]   28.45  28.20 1.90 1.20  26.48  31.03 1.01      837      233
##    m1[8]   31.36  30.97 2.48 1.51  28.97  34.55 1.01      427      217
##    m1[9]   34.26  33.75 3.10 1.85  31.32  38.19 1.01      319      205
##    m1[10]  37.16  36.53 3.73 2.18  33.67  42.04 1.01      271      182
##    y1[1]   11.04  11.46 4.46 3.52   4.29  16.96 1.01      629      371
##    y1[2]   13.86  14.13 4.43 3.39   6.95  20.00 1.01      762      282
##    y1[3]   16.86  17.10 3.81 3.09  10.82  22.36 1.01     1111      708
##    y1[4]   19.85  19.94 3.84 3.08  14.00  25.46 1.00     1259      856
##    y1[5]   22.64  22.60 3.76 3.08  17.20  28.63 1.00     1718     1324
##    y1[6]   25.59  25.44 4.12 3.10  19.93  31.55 1.00     2058     1242
##    y1[7]   28.52  28.26 4.22 3.27  22.76  34.84 1.01     1542      738
##    y1[8]   31.36  31.03 4.39 3.31  25.85  38.01 1.00     1282      398
##    y1[9]   34.18  33.69 5.18 3.49  28.10  40.97 1.00      828      344
##    y1[10]  37.01  36.35 5.31 3.64  30.96  44.13 1.01      597      390

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
data=list(N=n,
          x=sample(0:1,n,replace=T),
          y=sample(0:1,n,replace=T))

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -17.2496 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5       -13.393   0.000174194   1.81971e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -13.39
##      p        0.72
##      se       0.56
##      sp       0.64
##      ppv      0.79
##      npv      0.36
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -19.36 -19.00 1.39 1.17 -22.01 -17.83 1.00      776      876
##      p      0.49   0.48 0.29 0.37   0.05   0.95 1.01      693      566
##      se     0.54   0.54 0.14 0.15   0.30   0.78 1.00     2394     1653
##      sp     0.61   0.62 0.14 0.15   0.38   0.82 1.00     2947     1419
##      ppv    0.55   0.57 0.29 0.38   0.06   0.97 1.01      724      627
##      npv    0.56   0.59 0.29 0.37   0.06   0.97 1.01      726      548

2x2 cross table

Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)

n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)

p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n

Cramer'V  (chi2/n/(min(row,column)-1))^.5
  in 2x2  
  crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
      =(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5

kappa coefficient   k=(po-pe)/(1-pe)
  po: Observed agreement (proportion of times both raters agreed)
  pe: Expected agreement under independence
      po=p00+p11
        =(n0(1-p0)+n1p1)/n
      pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
        =n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/n

ex16-1.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
}

ex16-2.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
  real crv;
  crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
  real k;
  real po;
  po=(n0*(1-p0)+n1*p1)/n;
  real pe;
  pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
  k=(po-pe)/(1-pe);
}
n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -48.2981 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -34.2445    0.00167309   0.000743161           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -34.24
##      p0       0.17
##      p1       0.47
##      RR       2.80
##      OR       4.38
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -38.62 -38.31 1.06 0.74 -40.79 -37.61 1.00      946     1146
##      p0     0.19   0.18 0.07 0.07   0.09   0.31 1.00     1380      995
##      p1     0.47   0.47 0.09 0.09   0.33   0.61 1.00     2037     1292
##      RR     2.93   2.58 1.48 1.05   1.38   5.64 1.00     1423     1217
##      OR     4.93   4.07 3.38 2.27   1.65  11.35 1.00     1487     1191

data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -52.7725 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -34.2445   0.000159391   6.52948e-05      0.9447      0.9447        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -34.24
##      p0       0.17
##      p1       0.47
##      RR       2.80
##      OR       4.37
##      crv      0.32
##      k        0.63
##      po       0.65
##      pe       0.05
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -38.62 -38.31 1.06 0.74 -40.79 -37.61 1.00      946     1146
##      p0     0.19   0.18 0.07 0.07   0.09   0.31 1.00     1380      995
##      p1     0.47   0.47 0.09 0.09   0.33   0.61 1.00     2037     1292
##      RR     2.93   2.58 1.48 1.05   1.38   5.64 1.00     1423     1217
##      OR     4.93   4.07 3.38 2.27   1.65  11.35 1.00     1487     1191
##      crv    0.30   0.31 0.12 0.11   0.12   0.49 1.00     1615     1415
##      k      0.62   0.62 0.06 0.06   0.53   0.72 1.00     1671     1415
##      po     0.64   0.64 0.06 0.06   0.55   0.73 1.00     1700     1404
##      pe     0.05   0.05 0.00 0.00   0.05   0.06 1.00     1910     1637

point of subjective equality PSE

PSE: 50% threshold for sensing the difference between two stimuli is equal
JND: Just noticeable difference, difference between 50% threshold and 75%

r~B(n,p) #reaction for stimuli
p=1/(1+exp(-(a+b*x)))
x=x1-x0, x0,x1 is strength of stimuli

PSE=-a/b
JND=(log(0.75/0.25)-a)/b-PSE

ex6-3-0.stan

mulit logistic regression

data{
  int N;
  int m;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~binomial_logit(m,b0+b1*x);
}
n=30
m=10
x=runif(n,-2,2)
y=rbinom(n,m,1/(1+exp(-(-2+3*x))))

glm(y/m~x,family=binomial('logit'))
## 
## Call:  glm(formula = y/m ~ x, family = binomial("logit"))
## 
## Coefficients:
## (Intercept)            x  
##       -1.79         2.84  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  28 Residual
## Null Deviance:       22.5 
## Residual Deviance: 2.72  AIC: 14
data=list(N=n,m=m,x=x,y=y)

mdl=cmdstan_model('./ex6-3-0.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -694.414 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        8      -83.3571    0.00648888    0.00171875           1           1       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -83.36
##      b0      -1.79
##      b1       2.84
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -84.33 -84.05 0.95 0.68 -86.30 -83.42 1.01      784      850
##      b0    -1.84  -1.82 0.29 0.29  -2.36  -1.39 1.00      630      845
##      b1     2.92   2.89 0.36 0.36   2.34   3.57 1.01      612      870

b0=mcmc$draws('b0')
b1=mcmc$draws('b1')

pse=as.vector(-b0/b1)
quantile(pse,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)
##    0%  2.5%    5%   25%   50%   75%   95% 97.5%  100% 
## 0.426 0.508 0.522 0.587 0.630 0.676 0.741 0.762 0.861
jnd=(log(0.75/0.25)-b0)/b1-pse
quantile(jnd,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)
##    0%  2.5%    5%   25%   50%   75%   95% 97.5%  100% 
## 0.262 0.297 0.308 0.349 0.380 0.413 0.469 0.485 0.562